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Brownian dynamics simulations of bidisperse hard discs moving in two dimensions 
in a given steady and homogeneous shear flow are presented close to and above the 
glass transition density. The stationary structure functions and stresses of shear 
melted glass are compared quantitatively to parameter free numerical calculations 
for monodispcrsc hard discs using mode coupling theory (MCT) within the in- 
tegration through transients (ITT) framework. Theory qualitatively explains the 
properties of the yielding glass but quantitatively overestimates the shear driven 
stresses and structural anisotropies. 
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1. Introduction 

Shear flow can drive dense colloidal dispersions into states far from equilibrium. 
Especially of interest is the possibility to shear melt colloidal solids, in particular 
metastable colloidal glasses and gels, and to investigate shear-melted (yielding) 
colloidal glasses. Does a yield stress and/or yield strain exist (Petekidis et al. 2002)? 
Are shear-melted states necessarily heterogeneous (e.g. shear banded)? Does ageing 
prevent stationary states under steady shearing? Do hydrodynamic interactions 
dominate the properties in flow? These are just some of the questions whose answers 
will provide insights into the still murky glassy state. 

Recently the integration through transients (ITT) approach has been used to 
generalize the mode coupling theory (MCT) of the structural glass transition to the 
case of steady shearing (Fuchs Gates 2002, 2009), and to arbitrary homogeneous 
flows (Brader et al. 2008), albeit based on a number of approximations. First, 
homogeneous flow profiles are assumed from the outset, second, hydrodynamic in- 
teractions are neglected, and third, the mode coupling (mean-field-like) decoupling 
approximation, splitting a four point density fluctuation function into the prod- 
uct of two-point density correlators, has been applied. While simplified (schematic) 
models of MCT-ITT (Fuchs Gates 2003) can well be fitted to the hnear rheology 
(frequency dependent storage and loss moduli) and non-linear rheology (flow curves, 
viz. stress as function of shear rate) of model dispersions (Siebenbiirger et al. 2009), 
no fully quantitative comparison between the theory and data without additional 
approximations has been presented up to now. 
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We provide this first quantitative comparison, solving the MCT-ITT equations 
for hard discs confined to move in a plane, and performing Brownian dynamics 
simulations for a binary hard disc mixture, also in two dimensions. Computational 
(memory) constraints force us to work in two dimensions, which is still of some 
experimental relevance, as glasses have been observed in quasi-two-dimensional dis- 
persions (Bayer et al. 2007). Use of a binary mixture suppresses the nucleation rate 
sufficiently for even our long simulation runs to remain free of crystal nuclei. As 
hydrodynamic interactions are absent also in the simulation, and as Lees-Edwards 
boundary conditions (LE) in combination with the thermostat impose a homoge- 
neous shear rate, the integration of the equations of motion in the simulation ensure 
a stringent test of the theory. Results for the stationary shear and normal stresses, 
and for the shear-distorted microstructure are presented and discussed. The input 
quantities entering the theory, and transient density correlation functions, which 
(crucially) enter during intermediate MCT-ITT steps, are characterized also. 

Our work bears some similarity to the study by Miyazaki, Reichman, and Ya- 
mamoto, who, however, concentrated on fluid states under shear, and on their time 
dependent fluctuations (Miyazaki et al. 2004). We focus here on the yielding glass 
state, and its stationary, time-independent averages, which are not accessible to 
the theory in Ref. (Miyazaki et al. 2004). In order to compare quantitatively with 
theory, we aimed for better statistics than in comparable previous two-dimensional 
simulation studies of sheared glasses (Yamamoto Onuki 1998; Furukawa et al. 2009). 



The basic concept of the algorithm has been described in detail in three dimensions 
in (Scala et al. 2007) and can easily be adapted to two dimensions x and y. In 
order to prevent crystallisation at high densities we consider a binary mixture with 
a diameter ratio of D2/D1 = 1.4 at equal number concentrations (iVsmaii = -^big)- 
N ~ 1000 hard discs move in a simulation box of volume V with periodic Lees 
Edwards boundary conditions at packing fraction ip = ^ ■ [D\ + D\). The mass 
of the particles is set equal to unity, thermal energy fc^T sets the energy scale, and 
Di = D \s used as unit of length in the following. We choose our coordinate axes 
such that flow is in the x-direction and the shear gradient 7 is in the y-direction. 
After putting the particles on their initial positions we provide Gaussian distributed 
velocities which will be overlain by the linear shear flow. To propagate the system 
forwards in time we employ a semi-event-driven algorithm. For every particle at 
the time Iq the algorithm determines the possible collision time t with any other 
particle. This is easily achieved by solving the equation 



where Di and Dj denote the diameters of the particles i and j, the relative vector 
between both particles, and Vij the relative velocities. The smallest solution r,„i„ 
for all particle pairs determines the next event in the algorithm. All particles can 
then be propagated according to: r{Tmin) = '"(^0) + v T,nin- With the conservation 
of energy and momentum the binary collision laws impose new velocities v* for the 
particle j and i 



2. Simulation 



^{D, + D,) = \n 




(2.1) 



V* = Vj + {fij ■ Vij)r, 



V* =Vj - {f,j-Vij)f,j. 



(2.2) 
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Due to the boundary conditions any particle in the vicinity of the box-boundary 
can colhde with an image particle coming from the other end of the box. The boxes 
are constructed in such a way that they are translated with the velocity jL, where 
L denotes the size of the box. The boundary conditions thus have the consequence 
that the velocities tend to increase. Therefore a thermostat has to be introduced. 
After a time tb = 0.01 a so-called Brownian step sets in, which assures that the 
particles move diffusively for longer times. In the Brownian step at the time tb, all 
particle velocities are freshly drawn from a Gaussian distribution with {\v\'^) = 2. 
After that linear shear flow is imposed so that {v) = '^y{TB)i holds. 

As the system starts from a cubic lattice it is necessary to wait for the sys- 
tem to relax before meaningful stationary averages can be taken. The quantity of 
interest in the work presented here is the potential part of the stress: (Japii) — 
F^SiJ■(-^(O^J)a(''(Ou)0)^^^ with the relative force components (in a-direction) of 
particles i and j {F{t)ij)a, and the particles' relative distance components ii'{t)ij)0. 
The kinetic part will play no role, and thus has already been omitted. As we consider 
hard particles the forces must be calculated from the collision events. By observing 
the collisions within a certain time window Atc forces may be extracted using the 
change of momentum which occurs during the observation time. This leads to the 
evaluation algorithm (Lange et al. 2009; Foss Brady 2000) 

(Ta/3(7)=(^ J2 {^'"itch)a{r{tch)(3\ , (2.3) 

where the summation is over all collisions of particles i and j at time tc within the 
time window Atc. The procedure effectively sums the momentum changes {Avij)a 
in the a— direction multiplied by the relative distance of the particles {rij)p in the 
/?— direction. Here and below the brackets {■■■)s denote the average over different 
simulation runs. 

Additionally the shear stress can be computed via the contact value g{r): 

a^yii)^^^'- ^ fd9cos{0)sini0)g^^{r,0) (2.4) 

where g^^{r,6), g^^{r,6)... denote the 0-dependent partial contact values of the two 
components and dii,di2, ... the minimal distance between two particles. 6 is the 
polar angle. 

oo 

The Grcen-Kubo relation rj^y = ^ ipy / {<Jxy{t)crxy{t + s))ds holds for the non- 

^ 

sheared system. Thus the shear viscosity can be extracted from the simulation via 
(Alder et al. 1970) 




\ \coU J 



= 2fc^/i"ii \ / = y ( I > ^^^^^^ I ) ■ (2.5) 

where the sum runs over all collisions. The second pivot in this analysis is the equal- 
time structure factor and its deviation from the quiescent system. Exploiting that 
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^q{i) ~ {'k j cxp(iq(ri — fj)))^'^^ is a real quantity, we can extract it via 

^■jW^V^ E cos(q(r,-r,))\ (2.6) 

where the double sum runs over all pairs of particles i and j (NI) . The pairs ij are 
determined by the LE boundary conditions and the constraint of having the lowest 
distance among all possible image particles in the surrounding boxes. 
For low Peclet numbers (Peo=7£'^/-Do)in fluid states, g{r) can be expanded: g{r) = 
go{r)+2Pea^ gi{r)+0{PeQ) . This result can be used to derive the relative distortion 
of the structure factor in the linear response regime (Strating 1999). 

SS{q, 7 -> 0) = 2Peo J dr rgi (r) J dO cos{e) sin(6l)e[*(''«=" cos(e)+rg, sin(e))] ^3.7) 



while gi (r) can be obtained via 
27r 



dO g{r) cos{0) sm{0) ^ gi{r) + ©(Pc^ • xy/r^). (2.8) 



3. Mode Coupling Theory in the integrations through 
transients framework 

The MCT-ITT approach generalizes the MCT of the glass transition to colloidal 
dispersions under strong continuous shear. It considers N spherical particles with 
arbitrary interaction potential which move by Brownian motion relative to a given 
linear shear profile. An equation of motion for a transient density correlator $q(t) 
encodes structural rearrangements, and approximated generalized Green-Kubo laws 
relate stress relaxation to the decay of density fluctuations. 

The transient density correlator is defined by $q(i) ~ (jO*/9q(t) (i))/A^S'<j, where 

the density fluctuation is as usual Pq{t) ~ 12^=1 exp{iq- rj{t)}, and, in ITT, the 
average can be performed over the equilibrium Gibbs-Boltzmann ensemble. Thus 
the normalization of the initial value <I'q(0) = 1 is given by the equilibrium structure 
factor Sq. The time-dependent or shear- advected wavevector q{t) = {qx, Qy — '^tq^)^ 
appearing in the definition eliminates the affine particle motion with the flow field, 
and gives ^q{t) = 1 in the absence of Brownian motion. Shear flow coupled to 
random motion causes $q(i) to decay, as given by an (exact) equation of motion 
containing a retarded friction kernel which arises from the competition of particle 
caging and shear advection of fiuctuations 

$q(t) + rq(0 |$q(t) + ^ dt' mq{t, t' ) $q(t')| = , (3.1) 

where the initial decay rate contains Taylor dispersion as rq(t) = q^{t) / Sq^t)- The 
generalized friction kernel mq{t,t'), which is an autocorrelation function of fluctu- 
ating stresses, is approximated, following MCT precepts by an expression involving 
the structural rearrangements captured in the density correlators 

= Jwi^'^^mi^wF ^^^-^'^^^^-^''^ $k(.)(^-0$p(.)(i- 0(3.2) 
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with abbreviation p = q — k, n the particle density, and a vertex function given by 
V^qkp(t) = q(t) • {m Cfc(t) + p(t) Cp(t)) (3.3) 

where Ck is the Ornstein-Zernicke direct correlation function — {1 — 1/ Sk)/n. An 
additional memory kernel is neglected (Fuchs Gates 2009). The equilibrium struc- 
ture factor, Sk, encodes the particle interactions and introduces the experimental 
control parameters like density and temperature. Similarly, the potential part of the 
stress cTapii) = (fa^g )*•'''•' the non-equilibrium stationary state (neglecting the 
diagonal contribution that gives the pressure) is approximated assuming that stress 
relaxations can be computed from integrating the transient density correlations 



Flow also leads to the build up of shear-induced micro-structural changes, which, 
again integrating up the transient density correlators, can be found from 

S^ii) = S, + J^ dt <^l^_,^{t) . (3.5) 

A far smaller isotropic term in iS'q(7) (see (Fuchs Gates 2009)) is neglected here, as 
it is of importance for the plane perpendicular to the flow only. 



4. Numerical aspects 

The set of coupled Eqs. I3.1l3.2l and [3.3l was solved self-consistcntly using modifica- 
tions of standard algorithms and the Ng-iteration-scheme (Ng 2009). The functions 
were discretized on a 2D-Fourier-grid consisting of 101 grid points in either coor- 
dinate direction using a cutoff of qD = 30. More explicitly the Fourier-grid was 
discretized as qx,yD ~ —30, —29.4, . . . , 0, . . . , 29.4, 30 with D being the particle di- 
ameter. In order to enhance the accuracy, the advected direct correlation function 
Cq{t) in Eq. 13.31 was calculated from the input structure factor once it has been 
advected beyond the Fourier-grid cutoff up to qD = 100. For the time-integration 
of the convolution integral in Eq. 13.11 an initial time step of 10~^D^/r'o was used. 
The time-integration in the calculation of derived quantities Eqs. 13.41 and 13.51 was 
performed by dividing the integration interval into two sub-intervals, the first one 
containing times shorter, the other one times longer than t* = 1/(10007) ZJ^/Dq. 
The integration was then carried out on the intrinsic quasi-logarithmic grid of the 
MCT equation solver for times < t < t* and on an equally-spaced linear time-grid 
with spacing At = t* for times t > t* . 



5. Equilibrium and transient quantities 

The MCT-ITT approach uses structural equilibrium correlations as input, computes 
transient structural density correlators to encode the competition between flow- 
induced and Brownian motion, and calculates all stationary properties from time- 
integrals over the transient fluctuations. In this section, the input and intermediate 
quantities of the theoretical calculations are presented and discussed. 
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(a) Equilibrium structure factor 

The equilibrium structure factor Sq is the only input quantity to the MCT-ITT 
equations. It varies smoothly with density or temperature, but leads to the transi- 
tion from a shear-thinning fluid to a yielding glassy state at a glass transition density 
fc- Figure [T] shows modified hyper-netted chain structure factors of monodisperse 
hard discs va d — 2 from (Bayer et al. 2007) used in the MCT-ITT calculations. 
For comparison, the averaged structure factors obtained from the binary mixture 
simulation are shown also. In both cases, the value of the glass transition density 
ipc. is included in the curves in Fig.[TJ and only smooth changes in 5*^ are noticeable. 
The short range order at the average particle distance, as measured in the primary 
peak of Sq, increases with densification. Differences in the structure between the 
mono- and the bidisperse system become appreciable beyond the primary peak, and 
especially beyond the second peak in Sq. Also the height of the primary peak in Sq 
at Lpc differs, indicating that the averaged structure factor in the bidisperse system 
does not well characterize the local structure and caging in the simpler system, and 
that a comparison with a bidisperse MCT-ITT calculation should be performed. 
As this is outside the present numerical reach, we proceed by comparing results for 
the characteristic wavevectors denoted qi to 94 in Fig. [1] The angular dependence 
of the anisotropic structure will be explored along the special directions indicated 
in the insets. 




5 10 qD 15 20 25 5 10 qD 15 20 



Figure 1. Left panel: modified hyper-netted chain structure factor, used as input quan- 
tity to MCT, for separation parameters e = —10^^ (red), 10~^ (black) and 10~^ (blue). 
The critical packing fraction is <^c = 0.69899. The selected wavevectors are qiD = 3.0, 
q2D = 6.6, gs-D = 9.0 and q^D = 12.6. Right panel: averaged structure factors for the 
simulated binary system. The packing fractions are — 0.72, 0.74, 0.76 and ip = 0.79. 
The approximate critical packing fraction is Lpc — 0.79 ± 0.005, as estimated from the flow 
curves of Fig. [5l Both insets exhibit the magnitude and the relative orientation of the 
selected wavevectors and introduce the colour code, which is used later on throughout the 
comparison. The wavevectors for the right panel are: qiD = 2.8, g2-D — 5.4, gsD = 7.5 
and g4-D = 10.1. 



(b) Transient coherent density correlation functions 

The central quantities in MCT-ITT, which encode the competition between 
shear driven motion and random fluctuations, are the transient density correlators. 
Structural rearrangements manifest themselves as a (second) slow relaxation process 
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in the $q(i)j whose relaxation time depends sensitively on the distance to the glass 
transition, measured by the (relative) separation parameter e — [ip — ipc)/^c, and on 
the magnitude of the shear rate. Figured] shows representative curves in a fluid state 
(left column) and in a shear-melted glassy state (right column) for the wavevectors 
and directions defined in Fig. [TJ and for nine different shear rates. The denoted 
bare Peclet numbers measure the shear rate compared to the Brownian diffusion 
time estimated with the single particle diffusion coefhcient Dq at infinite dilution. 
Characteristically, for the present strongly viscoelastic system, shear affects the 
structural relaxation already at extremely small bare Peclet numbers Peo <SC 1 . The 
short time motion, which corresponds to the local diffusion of the particles within 
their neighbour cages, however, is not much affected as long as Pcq < f holds. 

In the fluid (e < 0), a linear response regime, where shear does not affect the 
decay of thermal fluctuations, is observed for small (dressed) Peclet or Weissenberg 
numbers Pe < 1, where Pe = jr measures the shear rate relative to the intrinsic 
relaxation time. The latter can be estimated from the relaxation of ^q2{t), viz. at 
the primary peak of Sq, where the structure relaxes most slowly in the quiescent 
case. The final (or a-) relaxation time t increases strongly when approaching the 
glass transition; with r oc (-e)^-^* predicted by MCT (Bayer et al. 2007). For the 
fluid state in Fig. [5J only Peo < f 0~^ is small enough that Pe < 1 holds and the 
final relaxation curves agree for different shear rates. In the glass (e > 0), shear is 
the origin of the decay of the otherwise frozen-in density fluctuations, and all nine 
shear rates lead to different final decays. While for small and larger wavevectors 
the final decay is rather isotropic, around the main peak in Sq some anisotropy is 
noticeable in the transient correlators. In the direction perpendicular to the flow, 

= (red), shear is not as efficient in decorrelating the density as in the other 
directions, which fall closer together. 

One of the central predictions of MCT-ITT concerns the existence of a scaling 
law describing the yielding of glassy states, which manifests itself by an approach to 
a master function ^q{t) ^q{t) for decreasing shear rate, 7 and e > 0, where 
the rescaled time t = 'jt agrees with the accumulated strain. This scaling law, which 
is quite obvious in Fig. [2l is tested quantitatively for some wavevectors in Fig. [3l 
and holds extremely well (for gi), or with pre- asymptotic corrections (for 52)- The 
shapes of the yielding master functions can be very well fitted with compressed 
exponentials, $q(i) — AqCxp (— (t/Tq)^i), but the parameters, including the ex- 
ponent /3q, depend on wavevector and orientation; see Tabled] for representative 
values. 

The quite good fit to compressed exponentials is at present only a numerical 
observation. Yet, close to the glass transition asymptotic expansions are possible 
which analytically predict the initial part of the yield master functions (Fuchs Gates 
2002,2003,2009): 



|.q(t ^ 0, £ = 0) ^ - hq ^C(7)/(A - -)) t + . . . . 

As the critical glass form factor and the critical amplitude hq are isotropic, 
this result suggests a rather isotropic yielding process right at the glass transition 
density. Approximating <I'q(t) by an exponential, it also provides an estimate for 
the final relaxation time under shear. Except for c^'^'> all quantities in the above 
formula have been determined in quiescent MCT (Bayer et al. 2007), and our values 
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10'* 10"' 10"- 10"' 10° lo' 10^ lo' 10* 10-'^ lo'' 10"* 10"^ 10° 10^ 10* lO'' 10* lO' 

t [D-/Do] t [D-/Do] 



Figure 2. Transient density correlators in MCT in the liquid state (left column) for 
e — —10""^ and in the glassy state (right column) at e = 10~^ at gi, 92,93 and 54 (from 
top to bottom) and Peo — 10~^ . . . 10~^ (from left to right). The colour code gives the 
orientations according to Fig. [1] The difference of the plateau heights in the second and 
fourth row emerges as a discrepancy of the magnitudes and orientations of wavevectors 
due to the selected discretization. 
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Figure 3. Top row: transient density correlators as functions of accumulated strain 'yt in 
the glassy state at e = 10"'^ for the wavevectors gi (top left), q2 (top right); collapse/ 
approach onto a master function appears for Peclet numbers ranging from Peo = 10"^ to 
10~^. Bottom right, master functions of this glassy state at gi (isotropic, lines) and q2 
(anisotropic, differently colored symbols and lines). The dashed black lines in this panel 
show compressed exponential fits $q(f) — Aq exp ( — (f/rq)^**) with parameters according 
to Table [1] Bottom left: rescaled final relaxation times Tq/j for = 0, Peo = 10~^ 
and £ = 10~^ 10"^ and 10"" (red from top to bottom) and lO"'^ (black). The blue line 

gives the estimate deduced from the stability analysis rq/7 ~ fq/{hq ^ c't"' /(A — |)) for 

— 3 5 'pjjg juget shows a close-up of the region around the first peak. 

(e.g. A = 0.698) differ only because of the coarser discretization in q space that is 
necessary under shear. Close to the glass transition, wc estimate c'^'^' « 3.5. As 
the bottom left panel in Fig. |3] and Fig. [2] both show, this estimate for the final 
relaxation time is quite good close to the transition, but deeper in the glass the 
correlators become somewhat more anisotropic. Specifically for qx ~ Q (red) with 
wavevector magnitudes around |q| ~ (72, the correlator slows down relative to other 
orientations. 

More detailed analytical predictions, are possible around where spatial and 
temporal dependences in the transient fiuctuations decouple (Fuchs Gates 2003): 

<f>q(t) - + K \fcW\ Q±{tlta, \i\ t.^c(7)/C|£|) + . . . . 

Q± is a universal function, which only depends on A, and t^j oc |e|~^'^^ is another 
MCT time scale that diverges at the glass transition (Bayer et al. 2007). This 
factorization property known from quiescent MCT generalizes to steady shear, and 
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q 


H/D 


q^/D 


qy/D 








qi red 


3.0 





3.0 


0.479 


1.121 


0.255 


qi green 


3.0 


2.4 


1.8 


0.481 


1.117 


0.263 


qi blue 


3.0 


3.0 





0.479 


1.112 


0.258 


qi magenta 


3.0 


1.8 


-2.4 


0.480 


1.117 


0.249 


q2 red 


6.6 





6.6 


0.856 


1.084 


0.811 


52 green 


6.5 


5.4 


3.6 


0.878 


1.384 


0.291 


52 blue 


6.6 


6.6 





0.875 


1.670 


0.363 


q2 magenta 


6.5 


3.6 


-5.4 


0.849 


1.617 


0.431 



Table 1. Parameters of ttie compressed exponentials j4q exp (— (t/rq)'^'') fitted to the yield 

master functions shown m Fig. [21 




t [D^/Dol t [D^/Dq] 

Figure 4. Factorization property: Xq(t) = ^'^(t'\^^%'^') ■ l^ft picture shows the liq- 
uid state at e = — 10^'',7 = 10~*, while the right picture depicts the glassy state at 
£ = 10^'', 7 = 10^**. Black curves give the solutions of the MCT equations for < qD < 24. 
At the times t' — 247 and t" = 1683 the correlators are approximately 0.01 above, re- 
spectively below their critical plateaus. Coloured curves are rescaled solutions of the mod- 
ified /3-scaling equation X{t) = g(y/tl\Zg\lJ/t}) with t^ = 965 and correspond to values 
g{7) _ I (green), 3.5 (red) and 5 (blue). The coefficient c'''' — 3.5 was also found in Fig. 
[3] and captures the shear driven decay of the glass structure. 



is an essential feature of the localization transition that underlies glass formation 
in MCT. In Fig. H the quantity Xq{t) = [*'„'') is plotted which should 

become wavevector and orientation independent if factorisation holds. This holds 
in the liquid, and in the glass, where however shear leads to strong (anisotropic) 
pre-asymptotic corrections already for 7^ > 0.01. 



6. Results and comparisons 

Based on the approximated generalized Green-Kubo relations of Eqs. p.4p and 
(|3.5p . and the properties of the transient correlators discussed in the previous sec- 
tion, MCT-ITT makes a number of predictions for stationary stresses and struc- 
tural correlations (Fuchs Gates 2002, 2009). In the following, we will test these 
qualitatively, but also quantitatively, by comparing them to the two-dimensional 
simulation data. 
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(a) Stationary stresses 

The quantity of most interest in nonlinear rheology is the shear stress <7xy{j). 
Flow curves giving the shear stress as function of the shear rate can be obtained in 
the simulations and are shown in Fig. [51 The viscosity decreases below its Newtonian 
value 7^0 = o'xyi'j upon increasing the shear rate already for small Fcq 

values. Shear thinning in which the stress increases less than linearly with 7, sets 
in at Pc of the order of unity, and this crossover shifts to lower and lower Pcq 
for increasing density. At the density around (pc « 0.79, the crossover leaves the 
accessible shear rate window. We use this as estimate of an ideal glass transition 
density, where the final relaxation time and the quiescent Newtonian viscosity {rjo oc 
r) diverge. Moreover, at this density the flow curves change from a characteristic 
S shape in the fluid, to exhibiting a rather 7-independent plateau at small shear 
rates. This change is the hallmark of the transition in MCT-ITT between a shear 
thinning fluid and a yielding glass (Siebenbiirger et al. 2009). 



9=0 J2 -'--+-'-'- 
9=0.74 -X 




Figure 5. Shear stress a^yij) versus shear rate given as Peclet number Peo = jD^/Do in 
MCT-ITT and simulation: The data points give simulation results for the densities denoted 
in the legend. Empty symbols show the results obtained via 1)2. 411 and filled symbols the 
results from (|2.3p . The Peo — >■ (black solid lines) were calculated according to (|2.5p . The 
solid lines show calculations in MCT-ITT for e = 10"^ 10"^ 10"", 10"^ -10"", -10"^ 
and -10"^ (from top to bottom; blue e > 0, red e ^ 0, green e < 0). MCT-ITT results 
are shifted downwards by a factor 0.1 to match the simulation results. 

The numerical MCT-ITT solutions show the same transition scenario, but be- 
cause of the approximations involved, exhibit a different critical density than the 
simulations. Even if theory and simulation considered the same (binary) system, 
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10"^ 10"^ 10"^ 10"^ 10"^ 10° 

Peo 

Figure 6. Viscosity 77(7) versus shear rate given as Peclet number Peo — -yD^/Do in 
MCT-ITT and simulation. The color codes are the same as in Figure [5] Again the empty 
symbols show the results of (|2.4p and filled symbols the results from H2.3I) while the 
Peo — >■ (black solid lines) were calculated according to H2.5|l . MCT-ITT results are 
shifted downwards by a factor 0.1 to match the simulation results. 

a difference in critical packing fraction would be expected as is well known also 
in three dimensions (Voigtmann et al. 2004). Thus the following MCT-ITT cal- 
culations match the relative separation e from tpc in order to compare with the 
simulations. As MCT-ITT is aimed at describing the long time structural motion, 
errors need to be anticipated in its description of short time properties. This is 
obvious in real dispersions, where hydrodynamic interactions (neglected in MCT- 
ITT) affect the short time diffusion coefficient Dg, but may arise in the following 
comparisons, also. A rescaling of the effective Peclet number Pcq*^ = ^D"^ jD^ would 
correct for this change in D^. In order not to introduce additional fit parameters we 
refrain from doing so, but anticipate that future comparisons may require ^ Dq. 

The approach to a yield scaling law, where the final decay of the transient 
correlators depends on the accumulated strain only (see Figs. [2] and [3]), predicts 
the existence of a (dynamical) yield stress cr^y(e) = <^xyici ~^ ^i^^^i which char- 
acterizes the shear melted glass. In the bidisperse hard disc mixture, at the glass 
transition it takes the (critical) value = cr^y(e = 0) « 0.3kBT/D^; below the 
glass transition, the yield stress jumps to zero, cr^^(£ < 0) = 0. Its quantitative 
prediction is quite a challenge for theory, because Eq. 13.41 shows that it requires 
an accurate calculation of the shear driven relaxation process. MCT-ITT overesti- 
mates the critical yield stress cr+jj° by roughly a factor 10 because, presumably, the 
decay of the transient correlators is too slow. Yet, of course the difference between 
the monodisperse system in the MCT-ITT calculation and the bidisperse simulated 
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system contributes in unknown way to the error. It appears reasonable to assume 
that mixing two species reduces the stresses under flow, which would explain part 
of the deviation. 

At larger shear rates, the flow curves from simulation appear to approach a second 
Newtonian plateau which, presumably, strongly depends on the hard-core character 
of excluded volume interactions and is outside the reach of the present MCT-ITT. 
The latter, by using its sole input Sq rather than pair potential, is not directly aware 
of hard-core constraints. We checked, however, that the states remain homogeneous 
and random up to the Peo values shown. 




10"' 10"' 10"* 10"-' 10"- 10"' 10° lo' 10"' 10"' 10"* 10"-' 10"- 10"' 10° lo' 



Pe„ Peo 

Figure 7. Left panel: Normal stress differences for the system at the densities 
given in the legend. The lines show MCT-ITT results at separation parameters 
e = 10^^,10"^,10~'',10~^-10~■*,-10"^ and -10"^ (from top to bottom) scaled by 
a factor 0.1 to match the simulation results. Right panel: Normal stress differences for the 
system divided by Peg. Color codes and scaling factors are the same as in the left panel. In 
both panels empty symbols show the results obtained via (|2.4p while filled symbols were 
calculated with (|2.3|l . 

Reassuringly, the same rescaling factor of 0.1 as for the shear stress brings 
theoretical and simulational normal stress differences, axx — c^yy, hrto register also; 
see Fig. [71 The normal stress differences are positive (indicating that the dispersion 
would swell after flowing through a nozzle), and show similar behavior to the stress, 
increasing like 7^ in the fluid for small shear rates, and leveling out onto a plateau 
in the yielding glass. 

(6) Distorted micro- structure 

The macroscopic stresses in the flowing dispersion are experimentally most im- 
portant, but provide only an averaged description of the local effects of shear. 
Spatially resolved information can be learned from the distorted structure factor 
6Sq{j) ~ Sqij) — Sq, which in MCT-ITT is connected to the stress via 

nksT f d'^q QiQj dcq ... 

as follows from Eqs. p.4p and p.5|) . 

Figure [5] shows colour-coded structure factors 5'q(7) as function of the two- 
dimensional wavevector q, with qx in the direction of flow and qy along the gradient 
direction. In the left column, panels with simulation data are compared to panels 
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Figure 8. Distorted structure factor ^(q, 7) from simulation (left column) and from theory 
(right column): for states in linear response regime in the liquid (top row: simulation 
ip = 0.78 and Pcq = 2 • lO"'* MCT e = -10"^ and Peo = W^*); in the glass at high shear 
rate (middle row: simulation ip = 0.79 and Peo = 2, MCT e = 10"^ and Peo = 10"^); 
and in the glass at small shear rate (bottom row: simulation (p = 0.79 and Peo = 2 ■ 10"'', 
MCT e = 10"^ and Peo = 10""). 

in the right column obtained in MCT-ITT. Scattering intensities are presented 
{i) (top row) in the linear response regime in the fluid, effectively measuring the 
equilibrium structure factors already shown in Fig. [TJ (ii) (middle row) at high 
shear in the glass, where all densities are in the shear thinning region, and (Hi) 
(bottom row) in the glass at low shear rate, where the yielding glassy state is 
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tested. While the fluid Sq{'j) is isotropic for small Peo (case (i)), as required by 
linear response theory, increasing Pcq to values around unity, leads to an ellipsoidal 
scattering ring, which is elongated along the so-called 'compressional axis' = 
—Qy, and more narrow along the 'extensional axis' = qy] this indicates that 
shear pushes particles together along the compressional and pulls particles apart 
along the extensional diagonals (Vermant Solomon 2005). Theory and simulation 
data in this representation qualitatively agree except for that MCT-ITT somewhat 
overestimates the anisotropic distortion of the glass structure at low Peo. 

A more careful look onto the distorted microstructure is provided by g-dependent 
cuts through S'q(7) along the directions colour-coded in Fig. [T] Important for the 
direction = (red) is the need, present in both simulation and theory, to av- 
erage >S'q(7) over a small but finite angle, because exactly at q^ — the structure 
oscillates noisily around zero. 

Especially of interest are the intensities of case (iii), where the stationary structure 
of the shear melted solid is studied. The density (p = 0.79 is not yet high enough to 
lie in the glass, but close enough to the glass transition so that the correlations at 
quite low rate, namely Peo = 2 10"'', closely resemble the ones of glassy states at 
very low shear rates. (While we can get good statistics for stresses at Peo = 2 10"^ 
for all densities, structure factors can not be sampled sufficiently there.) At this 
density oi ip = 0.79 also the equilibrium structure factor Sq can be obtained in 
long simulation runs, and the (relative) difference {Sq{'j,e) — Sq)/Sq can thus be 
determined. It is shown in Fig. IH] for two states in order to investigate the distorted 
structure of the shear melted glass in detail. We base the following discussion on 
the hypothesis that Peo = 2 • 10~^ a,t ip = 0.79 captures a glass-like state in the 
limit of low bare Peclet number. 

Figure [9] provides a sensitive test of the accuracy of the theoretical predictions. 
The lower left panel shows that the structure factor at vanishing shear rate Sq{'j — )■ 
0) jumps discontinuously at the glass transition; while 5'q(7 — >■ 0) — > 5'g holds in the 
fluid, 5'q(7 0) ^ Sq holds in the glass. Relative deviations (5'q(7 — 0) — Sq)/Sq of 
20% remain. Simulation finds quite isotropic deviations which show a maximum on 
the low-g side of the primary peak in 5"^. MCT-ITT predicts the absence of a linear 
response regime in >5'q(7) as function of the shear rate in the glass, and derives it 
from the existence of the yield scaling law in the transient correlators. Because 7 
sets the time scale for the final relaxation into the stationary state, the limit 7 ^> 
does not agree with the quiescent result 7 = 0. 

Quantitatively, MCT-ITT overestimates the distortion again by a factor up to 10, 
and finds a noticeable anisotropy, as discussed in context with Fig. |8l While the 
difference between the bidisperse and the monodisperse system may influence the 
comparison, we believe that the major origin of the error is that MCT-ITT under- 
estimates the speeding up of structural rearrangements caused by shear. The too 
slow transient correlators thus become anisotropic because the accumulated strain 
7 1 becomes too big before structural correlations have decayed. 

Qualitatively, aspects of the anisotropy predicted by MCT-ITT can be seen in 
the simulations at only slightly larger shear rates, like at Pcq = 0.2 shown in the 
middle left panel of Fig. [HI While along the two axis- and the extensional diagonal 
direction, the low-q wing of the primary peak in 5'q(7) becomes enhanced under 
shear, along the compressional axis (magenta) it gets suppressed, and the high-g 
wing is pushed up. Simulation also finds a suppression of the peak height along all 
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Figure 9. Direct comparison of the quantity (^(q, 7) — Sq)/Sq between simulation and 
MCT. The selected orientations are: = Q (red), — qy (green), qy = (blue) and 
qx — —qy (magenta) for tp = 0.79 in the simulation (left column). The Peo numbers are 
4 ■ 10"" (top row), 2 ■ 10"^ (middle row) and 2 ■ 10"^ (bottom row). For MCT (right 
column) the values are: £ = — 10~^, Peo = 10~* (top row), e = 10"'', Peo = 10~^ (middle 
row) and e — 1Q~^, Peo = 10~* (bottom row). Note that in the linear response case (top 
row) results were normalized by the shear rate 7. All structure factors are averaged over a 
small angle in the q direction (for an absolute |q|-value and an angle ip a small section Aq 
and Aip is defined; all q values lying within the defined section are used for the averaging. 

A(/9sim = 5°, A|q,hn| = 0.4,A<^MCT = 15°, A|qMCT| = 1.0 .). 



directions, which MCT-ITT reproduces along the diagonal directions. Overall the 
anisotropy and the magnitude of the distortions predicted by MCT-ITT remain too 
large, but the deviation decreases. The differences between the equilibrium structure 
factors Sq in the simulated and in the calculated system should be taken into 
account in future work, and presently preclude comparisons at larger wavevectors. 
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7. Conclusions and Outlook 

The first fully quantitative solutions of the MCT-ITT equations for the distorted 
niicrostructure and the stresses in steadily sheared two-dimensional shear-thinning 
fluids and yielding glasses of Brownian hard discs exhibit all the universal features 
discussed within schematic MCT-ITT models (Fuchs Gates 2003). They compare 
qualitatively well, but quantitatively with appreciable errors, with Brownian dy- 
namics simulations of a bidisperse mixture without hydrodynamic interactions in a 
linear shear profile, which for all states considered remains in an homogeneous and 
disordered state. The non-analytic behavior of the stationary properties, and the 
lack of a linear response regime throughout the (shear-melted) glass state, predicted 
by theory, can be found in the simulation. 
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